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The lattice Boltzmann method has been successfully applied for the simulation of flow through 
porous media in the creeping regime. Its technical properties, namely discretization, straightforward 
implementation and parallelization, are responsible for its popularity. However, flow through porous 
media is not restricted to near zero Reynolds numbers since inertial effects play a role in numerous 
natural and industrial processes. In this paper we investigate the capability of the lattice Boltzmann 
method to correctly describe flow in porous media at moderate Reynolds numbers. The selection 
of the lattice resolution, the collision kernel and the boundary conditions becomes increasingly 
important and the challenge is to keep artifacts due to compressibility effects at a minimum. The 
lattice Boltzmann results show an accurate quantitative agreement with Finite Element Method 
results and evidence the capability of the method to reproduce Darcy's law at low Reynolds numbers 
and Forchheimer's law at high Reynolds numbers. 

PACS numbers: 47.ll.-j 91.60.Np 47.56. +r 



I. INTRODUCTION 

Predicting the transport properties of porous media, 
like the fluid permeability, defined by Darcy's law, or 
heat conductivity, is of paramount importance in chemi- 
cal, mechanical, geological, environmental and petroleum 
industries. Flow situations in porous media are not 
restricted to the creeping flow regime, i.e., near zero 
Reynolds numbers where Darcy's law applies [1]. Many 
important natural and industrial processes are charac- 
terized by large Reynolds numbers, where inertial effects 
also play a role. Examples include gas flow through a 
catalytic converter, groundwater flow, filtration processes 
and the flow of air in our lungs [2] . Hence, accurate simu- 
lation methods arc needed to improve our understanding 
of these processes. 

The lattice Boltzmann (LB) method has become an 
efficient tool [3-5] as an alternative to a direct numerical 
solution of the Navier-Stokes equation [6, 7] for simulat- 
ing fluid flow in complex geometries such as porous me- 
dia. Historically, the LB method was developed from the 
lattice gas automata [5, 8], replacing the number of parti- 
cles in each lattice direction with the ensemble average of 
the single particle distribution function, and the discrete 
collision rule with a linear collision operator. In the LB 
method all computations involve local variables so that it 
can be parallelized easily [7] . Together with uniform grids 
and thus straightforward discretization, the LB method 
has become very popular in the field of porous media flow 
simulations. With the advent of more powerful comput- 
ers it became possible to perform detailed simulations of 
finely resolved complex samples [3, 4, 7, 9-15]. 

In this work we investigate the accuracy of the LB 



method in flow regimes beyond Darcy's regime. We focus 
on limitations of the method with respect to the lattice 
resolution and the selection of parameters. In particular, 
we address the requirement of reducing undesired com- 
pressibility effects by keeping the Mach number low, and 
how this influences the achievable maximum Reynolds 
numbers. In order to keep the compressibility as mod- 
erate as possible a new simulation setup combining an 
injection channel, density boundary conditions, and an 
external body force for driving the fluid is introduced. 
Nevertheless, for high Reynolds numbers, a low fluid vis- 
cosity and high resolution are required. Thus, a thorough 
investigation of the impact of the compressibility on the 
measured permeability is presented. 

To confirm the validity of the LB results, they are 
compared to Finite Element Method (FEM) simulations 
that are performed with the commercial software pack- 
age ANSYS. During the past decades, the FEM has 
been widely used to simulate fluid flow through porous 
media. It is known that FEM can deal with complex 
pore geometries and boundary conditions, see for ex- 
amples Refs. [16-18]. Tezduyar et al. [19] have devel- 
oped the so-called deforming-spatial-domain/space-time 
(DSD/ST) procedure for flow problems with deform- 
ing interfaces using the so-called Arbitrary Lagrangian- 
Eulcrian (ALE) methods and a space-time finite element 
method. This approach was based on fully resolved sim- 
ulations around particles and therefore computationally 
expensive in dense flows. For an overview of some finite 
element and finite difference techniques for incompress- 
ible fluid flow see Ref. [20] and for the efficiency of the 
solution algorithms Ref. [21]. In recent studies, Yazdchi 
et al. relate the macroscopic properties of porous media, 
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namely permeability and incrtial coefficients, to their mi- 
crostructure and porosity [22, 23]. 

The remainder of this article is organised as follows: In 
section II, an introduction to porous media flow includ- 
ing laminar and weakly inertial flow, i.e. creeping flow, 
named Darcy's flow, and Forchheimer's flow is given, re- 
spectively. In section III, the LB and FEM methods 
are introduced along with the description of the simu- 
lation setup. In section IV, we demonstrate quantitative 
agreement of LB and FEM simulations with theoretical 
predictions for fluid flow in porous media in the above 
mentioned regimes. Finally we analyze and compare the 
results before we conclude. 



where c is a positive dimcnsionlcss parameter and fh is 
the reference fluid density. The quadratic term in Eq. (3) 
describes a linear relationship between the solid-fluid hy- 
drodynamic interaction and the relative solid-fluid veloc- 
ity. Darcy's law is recovered in the limit of u s — >0. 

A transition from creeping flow to incrtial flow has 
been reported in many studies, see for example the work 
of Koch and Ladd on flow in a random array of cylin- 
ders [11] and spheres [12]. Also the existence of a tran- 
sition regime between the creeping and inertial regimes 
has been reported in the past [30]. This departure from 
Darcy's regime is of the order of u s 3 [22, 26, 31-33] and 
the relation between ((Vp)^) and u s in this short inter- 
val is then given by 



II. POROUS MEDIA FLOW 

Weak inertial flow, also named creeping flow (near zero 
Reynolds numbers) in porous media, can be described by 
Darcy's law [24, 25]. It is defined by 



((Vp) 2; ,)=-(^ + ^u s 2 | 

ft fi 



(4) 



where the cubic term in u s with positive parameter 7 
represents the weak inertia correction. 
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where the porous medium parameter k (always positive) 
is named the permeability of the medium. ((Vp) X2 ) rep- 
resents the average pressure gradient in the direction of 
the flow X2, see Fig. 1, and fi represents the dynamic vis- 
cosity of the fluid. u s represents the superficial velocity 
defined by 
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u 2 dV, 
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where 1*2 is the fluid velocity in the flow direction. Vt 
and Vf represent the total volume and the fluid volume, 
respectively. The average velocity within the porous 
medium (u) is related to the superficial velocity by u s = 
e(u), where e is the porosity of the medium. According 
to Eq. (1), Darcy's law corresponds to a linear relation 
between the average pressure gradient ((Vp)^) an d the 
superficial velocity u s , in the literature also referred to as 
seepage velocity [26]. Darcy's law was obtained empiri- 
cally in 1846, but it can be derived from the continuous 
moment and mass balance assuming that the solid-fluid 
hydrodynamic interaction is proportional to the relative 
solid-fluid velocity [27]. 

When the Reynolds number is increased, inertia ef- 
fects become relevant (Re ps 0(1)). The average pres- 
sure gradient {(Vp) X2 ) and the superficial velocity u s do 
not follow a linear relation anymore as it was empirically 
shown by Forchhcimcr [28, 29]. For this flow regime, 
named Forchheimer's regime, a quadratic term in u s is 
included to take into account the inertia effects. Accord- 
ing to Massarani [27], Forchheimer's law can be written 
as 
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III. SIMULATION METHODS 

A. The lattice Boltzmann Method 

The integration of the Boltzmann equation on a regu- 
lar lattice using a discrete set of velocities Cj defines the 
lattice Boltmann method. The lattice is defined by the 
spacing Ax, with the discrete velocity in units of Ax/ At, 
where At represents the timestep. Thus, the basic differ- 
ential equation for the method is 

/ i (x+Atc i ,t+At)-/ i (x,t) = -— (/i(x,t) - /rCM)) , 

T 

(5) 

where /i(x, t) represents the number of particles mov- 
ing at position x with velocity Cj at discrete time t. 
The term on the right hand side represents the colli- 
sion operator as introduced by Bhatnagar, Gross and 
Krook (BGK) [34, 35], which approximates the collision 
through a linearization towards the equilibrium distribu- 
tion function /° q (x, t) with a unique relaxation time r. 
This value is restricted to rj At > 1/2 assuring positive 
viscosity. If t/ At approaches 1/2 numerical instabilities 
can arise [36]. This collision operator is often referred 
to as "single relaxation time" or BGK model. Under 
the assumption of very low Knudsen and Mach numbers 
(Ma) which assures small compressibility effects, / { eq (x, t) 
is calculated by a second order Taylor expansion of the 
Maxwell distribution as [5] 



(u-Cj) S 

2c s 4 



2c s 2 



(6) 



p o is a reference density. The coefficients Wj are called lat- 
tice weights and are chosen to assure conservation of mass 
and momentum. They differ with lattice type, number of 
space dimensions and number of discrete velocities. For 
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the cubic fashion 3D lattice with 19 discrete velocities 
(i = 0, ..,18) used in this work they are 1/3, 1/18 and 
1/36 for the rest particles, the particles moving parallel 
to the xi, X2, or X3, and the particles moving in diagonal 
directions, respectively [5, 37]. Even though the system 
of interest in this paper is intrinsically two-dimensional, 
we apply a three-dimensional implementation of the LB 
method and compute the flow in a very flat simulation 
domain with periodic boundary conditions in the £3 di- 
rection. We do not expect this to have any influence on 
the simulation results, but it allows us to use our well 
tested implementation "LB3D". The only disadvantage 
are higher computational costs, but we do not report on 
the amount of CPU time required for a given flow prob- 
lem in this paper. 

No-slip boundary conditions on the solid walls are im- 
plemented by mid-plane bounce back rules [38]. The 
BGK model is known to suffer from an artificial viscosity 
dependent slip at boundaries if these boundary condi- 
tions are used. An alternative approach for the colli- 
sion operator, which reduces this well known drawback 
of the BGK model, is the multi relaxation time (MRT) 
method [13, 14]. Here, the right hand side of Eq. (5) is 
replaced by the expression 

-At [M- 1 • S • (m(x, t) - m oc »(x, *))] . , (7) 

where M is a linear transformation chosen such that the 
moments 

m i (x > t)=53M ij / j (x,*) (8) 
3 

represent hydrodynamic modes of the problem. We use 
the definitions given in [39], where mo is the fluid den- 
sity, m,2 represents the kinetic energy, m, with i = 3, 5, 7 
is the momentum flux and m,, with i = 9, 11, 13, 14, 15 
are components of the symmetric traceless viscous stress 
tensor. During the collision step the density and the 
momentum flux are conserved so that m; = with 
i = 0,3, 5, 7. The non-conserved equilibrium moments 
m° q , i 0,3, 5, 7, are assumed to be functions of these 
conserved moments and explicitly given in [39]. S is a 
diagonal matrix Sij = SiSij. The diagonal elements 
T; = 1 /§i in the collision matrix are the relaxation time of 
the moment m- % . One has §0 = S3 = §5 = §7 = 0, because 
the corresponding moments are conserved, si = l/rbuik 
describes the relaxation of the energy and sg = sn = 
S13 = S14 = S15 = 1/t the relaxation of the stress tensor 
components. The remaining diagonal elements of S are 
chosen such that one has 

S = diag(0, l/r bulk , 1.4, 0, 1.2, 0, 1.2, 0, 1.2, 1/t, 

1.4, 1/t, 1.4, 1/r, 1/r, 1/t, 1.98, 1.98, 1.98), { ' 

to optimize the algorithm performance [39, 40]. The two 
relaxation times r and Tbuik, are restricted as well as for 
the BGK method to be > At/2, but allow to define the 
kinematic and bulk viscosity, respectively. This multi 



relaxation time scheme is commonly referred to as "two 
relaxation time" (TRT) method. An alternative TRT 
implementation can be found in [41, 42]. 

The macroscopic density p(x, t) and velocity u(x, t) are 
obtained from /j(x, t) as 

p(x,i)= A ^/ i ( X ,i), (10) 

i 

u(x,i) = -A^^/ i (x,i)c i . (11) 

The pressure is given by 

p(x,t) =c s 2 p( X ,t), (12) 

where c s = l/\/3(Ax/At) is the speed of sound [5, 8]. 
The kinematic viscosity of the fluid v — /i/p is a function 
of the discretization parameters, Aa; and At, and the 
relaxation time t [43, 44]. It is given by 

c s 2 At 1 t \ , . 

^^2-{ 2 At- 1 )- ^ 

B. The Finite Element Method 

The velocity and pressure profiles through the system 
can be obtained from the solution of the conservation 
laws, namely, the continuity equation (conservation of 
mass) and the Navier-Stokes equations (conservation of 
momentum). In the absence of body forces, but assuming 
a constant density (i.e. incompressible flow) and steady 
state flow conditions, the equations of conservation of 
mass and momentum for a Newtonian fluid are simplified 
to 

V-u = 0, (14) 

p(u • Vu) = -Vp + /iV 2 u. (15) 

The FEM makes use of the variational formulations that 
allow the transformation of the above equations into 
a system of linear algebraic equations, which can be 
solved using a simple LU decomposition or iterative algo- 
rithms [45, 46]. Stable discretizations of the above equa- 
tions are difficult to construct and it is known that the in- 
comprcssibility constraint is not strongly enforced when 
using continuous polynomial shape functions for pres- 
sure. Sec [21] for a detailed theory and discussion. Lang- 
tangen et al. present an overview of the most common 
numerical solution strategies, including fully implicit for- 
mulations, artificial compressibility methods, penalty for- 
mulations and operator splitting methods [47] . Using the 
conventional FEM scheme, we solve the above equations 
with the commercial software ANSYS [48]. The nonlin- 
ear solution procedure used in ANSYS belongs to a gen- 
eral class of Semi-Implicit Methods for Pressure Linked 
Equations (SIMPLE). On the flow domain, the steady 
state Navier-Stokes equations combined with the conti- 
nuity equations are discretized into linear triangular ele- 
ments. They are then solved using a segregated sequen- 
tial solution algorithm. This means that element ma- 
trices are formed, assembled and the resulting system is 
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solved using the Gaussian elimination algorithm for each 
degree of freedom separately. The number of iterations 
required to achieve a converged solution may vary con- 
siderably, depending on the number of elements, inertial 
contribution and the stability of the problem. Some more 
technical details arc given in [48]. 

By knowing the fluid velocity field, the superficial ve- 
locity is then calculated from Eq. (2). Recently, using 
FEM simulations, Yazdchi and Luding [22] show that for 
both ordered and random fibre arrays, the weak inertia 
correction to the linear Darcy relation is third power in 
superficial velocity, up to small Re ps 1-5. When attempt- 
ing to fit the data with a particularly simple relation, a 
non-integer power law performs astonishingly well up to 
the moderate Re « 30. 

A typical unstructured, fine and triangular FEM mesh 
is shown in Fig. 1. The mesh size effect is examined by 
comparing the simulation results for different resolutions. 
The mesh refinement is done on element level, meaning 
that a finer grid is overlaid on the coarse one. In order 
to be able to apply periodic boundary conditions in x\ 
direction, see Fig. 1, we discretize the system such that 
we come up with the same number of nodes on the left 
and right boundary. Periodic boundary conditions are 
applied by setting additional constrains, i.e. same veloc- 
ity, on these nodes. 



C. Computational Domain and Simulation Setup 

The computational domain is presented in Fig. 1. The 
porous sample has a length of X$ representing 0.9 A2, 
and it is composed of 266 circles, randomly distributed 
with a radius r of 1.25% of X2. A minimum distance 
A m ; n of 3.125% of X2 is imposed between the circle cen- 
ters. The sample porosity then follows to be e — 0.63730. 

For the LB method three different resolutions are used, 
namely low resolution (L), intermediate resolution (M), 
and high resolution (H), see Table. I for details. For the 
FEM method also three resolutions are used, correspond- 
ing to 22048, 49670, and 982376 triangular elements, re- 
spectively. These meshes are referred to using the same 
abbreviations as for the LB simulations. However, one 
should keep in mind that the number of discretisation 
elements used in both methods cannot be compared eas- 
ily since the LB simulations utilize a regular Cartesian 
lattice, while the FEM simulations are based on unstruc- 
tured grids with locally varying resolution. 

In the LB simulations a constant acceleration g = ge^ 
drives the fluid using Guo's method [49]. This accelera- 
tion defines the pressure gradient and acts on the fluid 
inside the sample, i.e. from X2 = Xq to X2 = Xq + Xg. 
Together with this, on-site pressure boundary conditions 
are implemented on the inlet-plane X2 — and on the 
outlet-plane X2 = X2 setting a constant density value p 
(reference density, i.e. fi = 1) 011 both planes. Sec 
Eq. (12) for the relation between the density and pres- 
sure within the LB method. These pressure boundary 



conditions arc implemented using the method of Zou & 
He [50, 51]. Different values of the acceleration g are used 
in order to study different Re regimes. In our earlier work 
we proposed to use just on-site pressure boundary con- 
ditions at the inlet- and outlet-planes, (J] fi = 1 + 6 and 

fi = 1 — 8, respectively) to impose a pressure drop 
of 2(5 [15]. The here presented new setup allows to re- 
duce compressibility effects inside the sample when the 
Reynolds number is increased. Periodic boundary condi- 
tions are applied in the x\ direction. 

In the FEM simulations, a pressure drop is imposed 
while the density is kept constant to drive the fluid. Fur- 
ther, we impose zero velocity on the surface of the fibres 
and also apply periodic boundary conditions in the X\ 
direction. 



IV. SIMULATION RESULTS 

We present a detailed analysis of the simulation of 
porous media flow at low to moderate Reynolds num- 
bers using the LB method. Special attention is given 
to its collision kernel, the discretization and the value of 
the relaxation time. To increase Re either the average 
flow velocity can be increased or the fluid viscosity can 
be decreased. However, the range of these parameters 
is limited: the lattice Boltzmann method is known to 
reproduce the Navier-Stokes equations in the low Mach 
number limit only, i.e., at high flow velocities compress- 
ibility artifacts can occur and render the results invalid. 
To validate the obtained data and to understand the im- 
pact of these limits on the precision of the LB results, we 
present a quantitative comparison with FEM results and 
theoretical predictions. 

For the LB simulations 100,000 timesteps assure the 
steady state. The superficial velocity defined by Eq. (2) 
is calculated from the steady state LB data by 

where X represents all fluid nodes in the simulation do- 
main. The average pressure gradient is expected to be 

((Vp)x a )=-5ft. (17) 
The Mach and Reynolds numbers are defined by 

Ma=M (i 8 ) 

Rc=Ml (19) 

Fig. 2 shows the relative maximum fluid density as a 
measure for the compressibility. The data are obtained 
from LB-MRT simulations of flow in the model geome- 
try introduced above. While the open symbols represent 
data obtained using standard pressure drop boundary 
conditions (see [15]), the closed symbols describe data 
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FIG. 1. Computational domain with dimensions X\ x X2 (Xi = 0.4 Xi). The sample is composed of 266 circles with radius 
r representing 1.25% of X2. They are randomly distributed assuring a minimum circle center distance A m i n of 3.125% of Xi. 
Pure fluids chambers with length of 5% of Xi are placed after and before the sample. The right inset shows a zoom of a typical 
unstructured, fine and triangular FEM mesh. 



TABLE I. Domain discretization for the LB simulations. 



Resolution 


X\ Xi Xs Xc r A m i n 
Aa: Aa: Ax Ax Ax Ax 


low (L) 
medium (M) 
high (H) 


128 320 288 16 4 10 
256 640 576 32 8 20 
512 1280 1152 64 16 40 



obtained with the newly proposed simulation setup com- 
bining pressure boundary conditions with a body force 
g. In Fig. 2 the compressibility and thus the increase of 
the relative maximum density in the system increases for 
higher Re. Further, a clear reduction of compressibility 
effects is obtained when the newly proposed simulation 
setup is used. 

Fig. 2 also shows the effect of a reduced relaxation 
time (and thus viscosity) on the maximum reachable Re 
and the related compressibility effects. The compressibil- 
ity of the fluid starts to become less important for lower 
viscosities and the combination of a low relaxation time 
(r/At = 0.6) together with the combined pressure and 
body force boundary conditions leads to higher maximum 
Reynolds numbers. With this combination, Re of the or- 
der of 30 can be reached before, at large Re, fluctuations 
become too large for reliable analysis of the simulations. 

Fig. 3 shows the permeability estimate obtained from 
LB simulations using different discretizations and colli- 
sion kernels. In our previous work it was shown for BGK 
simulations of 3D Poiscuillc flow that when a relaxation 
time close to r/At = 0.8 is used, the dependency of the 
results on the discretization is minimal [14]. As it is also 
demonstrated in the same article, finding such an optimal 
value for the relaxation time is impossible for more re- 
alistic samples such as a discretized Fontainebleau sand- 



stone or the array of cylinders as used here: a strong de- 
pendency on the discretization appears when the BGK 
model is used, even though the relaxation time is set to 
t I At = 0.8. By using the MRT model this dependency 
can be decreased substantially. As can be observed in 
Fig. 3, the results of both collision kernels are qualita- 
tively similar, i.e., both methods correctly reproduce the 
expected shape of the curve for comparable maximum 
Re. One can see that in the case of the low resolu- 
tion sample (L) the departure from the Darcy's regime 
(constant permeability) is to higher permeability values, 
which is unphysical since c in Eq. (3) is a non-negative 
parameter. In the case of intermediate (M) and high res- 
olution (H), the departure is to lower permeability values, 
but only using the high resolution sample it is possible to 
simulate high-enough Reynolds numbers to analyze this 
phenomenon. 

Fig. 4 shows the permeability estimation for the 
LB-MRT method using different values for the relaxation 
time. Due the resolution used for this study (H) together 
with the MRT collision kernel it is not surprising that 
the results do not differ much quantitatively, but only by 
about 3%. The main difference is in the highest possible 
Reynolds numbers which can be reached by using a small 
value of the relaxation time r/At = 0.6, in accordance 
with Fig. 2. As mentioned above numerical instability 
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FIG. 2. Normalized maximum fluid density as a measure 
for the effect of compressibility versus Reynolds number. 
Open symbols denote data obtained from LB-MRT simula- 
tions where the fluid is driven by an imposed pressure drop 
(p-BC) as suggested in [15]. When combining a fixed density 
at the in- and outlet of the domain with a driving body force 
(p-BC + g), compressibility effects can be reduced substan- 
tially. Together with low values of the relaxation time r/At 
higher Re can be reached (solid symbols). 
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FIG. 3. Resolution and collision kernel (r/At = 0.6) anal- 
ysis for permeability estimation vs. Re. The results show a 
dependency with the domain discretization, which is stronger 
when the BGK collision kernel is used. 

arises when the relaxation time approaches t/ At = 0.5. 
For this reason the smallest value used in this work is 
r/At = 0.6. 

For porous media flow calculations using LB-MRT and 
t I At = 0.6, Fig. 2 shows that at high Reynolds numbers 
the compressibility is of the order of 30%. Furthermore 
as we can see in Fig. 5 the superficial velocity remains low 
enough to keep Mach numbers below rj 10 . However, 
inside the sample there are zones with high velocity. In- 
deed, for the two higher Reynolds numbers simulated, on 
~ 25% and ~ 45% of the fluid nodes the velocity is higher 




Re 

FIG. 4. Relaxation time analysis for permeability estimation 
vs. Re calculated using the LB-MRT method and the high 
resolution (H) computational domain. Different values for 
the relaxation time r/At were used. The results only differ 
by ~ 3% and the highest Re can be reached with the smallest 
value of r/At = 0.6. 

than 20% of the speed of sound c s , see Fig. 5. Such high 
velocities and strong compressibility which are present 
at high Reynolds numbers highly question the validity of 
the results. 



1.0 ^r-^ ^r-^ 1 n 1 — — 




FIG. 5. Superficial velocity u s and maximal velocity it™ ax 
inside the sample for LB-MRT and r/At = 0.6. The inset 
shows the percentage of fluid nodes with velocity u y higher 
than 20% of the speed of sound. 

To investigate the validity of the results FEM simula- 
tions are performed as a benchmark for the LB simula- 
tions. In Fig. 6 the data obtained using both methods 
are plotted. One can see that the FEM results also show 
a dependence on the sample resolution. As stated above, 
in the case of the FEM the L, M, and H consist of do- 
mains with 22048, 49670, and 982376 elements. The res- 
olution dependency is also demonstrated in the inset of 
the figure, where the permeability is shown as a function 
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of resolution for a Reynolds number of Re w 10 -3 . For 
simulations using 1 x 10 6 elements, there is a difference 
of fa 2% between LB and FEM results. This can be ex- 
plained by several factors: neither FEM nor LB results 
are fully converged with respect to the required number 
of elements, but performing a large number of simula- 
tions with substantially higher resolution is not feasible 
with the computational resources available - in particu- 
lar since the final values are expected to change by not 
more than a few percent. Since the LB results are sys- 
tematically lower than the FEM results, a possible expla- 
nation can be the loss of accuracy of LB due to the rel- 
atively small choice for the relaxation time (r/At = 0.6) 
as demonstrated in Fig. 4. In any case, the deviation is 
small and thus it can be concluded that the data in Fig. 6 
shows a qualitative and quantitative agreement when the 
Reynolds number increases. Fig. 6 shows that for both 
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FIG. 6. Comparison of the permeability vs. Re obtained using 
the LB-MRT (r/At = 0.6 and H) and different resolutions 
of the FEM grid. The inset shows the resolution dependent 
permeability for LB (LB-MRT and r/At = 0.6) and FEM 
simulations at Re w 10 -3 . 



methods the value of the permeability, beyond validity 
of Darcy's law, Eq. (1), drops consistently. Darcy's law 
is limited to the weak inertia regime, i.e. low Reynolds 
numbers. To analyze the validity of the simulated data 
for high Reynolds numbers, we plot u s vs. ((Vp) X2 ) m 
Fig. 7 using the H data and r/At = 0.6 for LB-MRT. 
For both methods, at small velocity values (small Re), 
the flow follows Darcy's prediction (constant slope) and 
Eq. (1) accurately fits the data. When the velocity in- 
creases the measured permeability departs from Darcy's 
law and the best possible fit is obtained using Eq. (4), 
which confirms the existence of a transition regime with 
cubic velocity dependence. Only for higher velocities the 
flow enters Forchheimer's regime and the data can be 
fitted accurately by Eq. (3). The departure from the dif- 
ferent flow regimes is clearly seen in the bottom panel 
of Fig. 7, where the relative error of the fits is plotted 



e{w) 



(20) 



where the indices s i m and fit represent the results ob- 
tained by simulation and from the fit, respectively. It can 
clearly be seen from both datasets that Darcy's regime 
holds until Re « 10° and the inertia transition correc- 
tion expression holds until Re « 10 1 . The Forchheimer's 
regime describes the flow for higher Re. 
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FIG. 7. (Top) Fit of the expressions which relate the super- 
ficial velocity u s and the average pressure gradient {(Vp) X2 ) 
for the porous media flow. (Bottom) Relative error of the fit 
presented above. The relative error shows that Darcy's law 
fits the simulation results for Re < 10° . Above this value the 
fit of the expressions given by Eq. (4) is better than the fit 
with Darcy's law. 



V. CONCLUSIONS 

We presented an analysis of the calculation of the flow 
in porous media in different flow regimes using the LB 
method and compared our results to FEM simulations. 



8 



Special attention has been given to discretization, selec- 
tion of a collision kernel and choice of parameters. For 
the LB method a simulation setup combining a constant 
pressure at the inlet and outlet and an external accelera- 
tion acting inside the sample has been proposed in order 
to reduce compressibility artifacts at high Re. However, 
at higher Re a substantial number of lattice nodes shows 
flow velocities beyond 20% of the speed of sound and 
compressibility effects can clearly be observed. In or- 
der to clarify the validity of these results we compare 
our data to FEM simulations and demonstrate a good 
quantitative agreement for the full range of Re studied. 
Furthermore, the data show good agreement with theo- 
retical predictions, demonstrating that the range of Re 
studied in this work is well accessible for the LB method 
and that compressibility effects only have a minor influ- 
ence. The results accurately predict three different flow 
regimes. These are Darcy's regime for Re < 10°, where 
the average pressure gradient ((Vp) X2 ) and the surface 



velocity u s obey a linear relationship. Secondly, a tran- 
sient regime 10° < Re < 10 1 , where the flow is mod- 
eled by Darcy's law plus an inertia correction term rep- 
resented by a cubic term on u s , see Eq. (4). Finally, for 
higher Reynolds numbers (Re > 10 1 ) the porous media 
flow follows Forchheimer's prediction up to the limit of 
both simulation methods of Re « 30. 



VI. ACKNOWLEDGMENTS 

Financial support by STW (STW-MuST program, 
project number 10120), NWO/STW (VIDI project 10787 
and VICI project 10828) and the FOM-Shell IPOG-II 
IPP is highly acknowledged. Computing time was pro- 
vided by the Julich Supcrcomputing Centre (JSC) and 
the Scientific Supercomputing Center Karlsruhe (SSCK). 



[1] J. Bear. Dynamics of fluids in porous media. Elsevier 

(New York), 1972. 
[2] K. Moutsopoulos, I. Papaspyros, and V. Tsihrintzis. 

Experimental investigation of inertial flow processes in 

porous media. Journal of Hydrology, 374(3-4) :242-254, 

2009. 

[3] B. Ferreol and D. Rothman. Lattice-Boltzmann simula- 
tions of flow through Fontainebleau sandstone. Transp. 
Porous Media, 20(1-2) :3-20, 1995. 

[4] A. Cancelliere, C. Chang, E. Foti, D. H. Rothman, and 
S. Succi. The permeability of a random medium: Com- 
parison of simulation with theory. Phys. Fluids., 2:2085- 
2088, 1990. 

[5] Y. H. Qian, D. d'Humieres, and P. Lallemand. Lattice 
BGK models for Navier-Stokes equation. Europhys. Lett., 
17(6):479-484, 1992. 

[6] N. Martys and E. Garboczi. Length scales relating the 
fluid permeability and electrical conductivity in random 
two-dimensional model porous media. Phys. Rev. B, 
46(10):6080-6090, 1992. 

[7] C. Manwart, U. Aaltosalmi, A. Koponen, R. Hilfer, and 
J. Timonen. Lattice-Boltzmann and finite-difference sim- 
ulations for the permeability of three-dimensional porous 
media. Phys. Rev. E, 66(1):016702, 2002. 

[8] S. Succi. The lattice Boltzmann equation for fluid dy- 
namics and beyond. Oxford University Press, 2001. 

[9] E. Aharonov and D. H. Rothman. Non-Newtonian flow 
(through porous media): a lattice Boltzmann method. 
Geophys. Research Lett., 20:679, 1993. 
[10] F. Auzerais, J. Dunsmuir, B. Ferreol, N. Martys, J. Ol- 
son, T. Ramakrishnan, D. Rothman, and L. Schwartz. 
Transport in sandstone: a study based on three di- 
mensional microtomography. Geophys. Research Lett., 
23(7):705, 1996. 
[11] D. Koch and A. Ladd. Moderate Reynolds number flows 
through periodic and random arrays of aligned cylinders. 
J. Fluid Mech., 349:31-66, 1997. 
[12] R. Hill, D. Koch, and A. Ladd. The first effects of fluid 
inertia on flows in ordered and random arrays of spheres. 



J. Fluid Mech., 448:213-241, 2001. 

[13] C. Pan, L.-S. Luo, and C. T. Miller. An evaluation of 
lattice Boltzmann schemes for porous medium flow sim- 
ulation. Computers & Fluids, 35:898-909, 2006. 

[14] A. Narvaez, T. Zauner, F. Raischel, R. Hilfer, and 
J. Harting. Quantitative analysis of numerical esti- 
mates for the permeability of porous media from lattice- 
Boltzmann simulations. J. Stat. Mech.: Theor. Exp., 
2010P11026, 2010. 

[15] A. Narvaez and J. Harting. Evaluation of pressure 
boundary conditions for permeability calculations using 
the lattice-Boltzmann method. Adv. in Appl. Math, and 
Mech., 2:685-700, 2010. 

[16] N. Zabaras and D. Samanta. A stabilized volume- 
averaging finite element method for flow in porous and 
binary alloy solidification processes. Int. J. for Num. 
Meth. in Eng., 60(5): 1-38, 2004. 

[17] V. Girault and P. Raviart. Finite elements approxima- 
tion of the Navier-Stokes equations. Springer Series SCM, 
1986. 

[18] F. Thomasset. Implementation of finite element methods 
for Navier-Stokes equations. Springer- Verlag, New York, 
1981. 

[19] T. Tezduyar, M. Behr, and J. Liou. A new strategy for 
finite element computations involving moving boundaries 
and interfaces-the deforming spatial-domain/space-time 
procedure: I. the concept and the preliminary numeri- 
cal tests. Comput. Meth. Appl. Mech. Eng., 94:339-351, 
1992. 

[20] C. Fletcher. Computational Techniques for Fluid Dynam- 
ics. Springer, 1991. 

[21] S. Turek. Efficient Solvers for Incompressible Flow 
Problems: An Algorithmic and Computational Approach. 
Springer, 1999. 

[22] K. Yazdchi and S. Luding. Towards unified drag laws for 
inertial flow through fibrous materials. Chemical Engi- 
neering Journal, 207:35-48, 2012. 

[23] K. Yazdchi, S. Srivastava and S. Luding. Micro-macro 
relations for flow through random arrays of cylinders. 



9 



Composites Part A, 43:2007-2020, 2012. 

H. Darcy. Les fontaines publiques de la ville de Dijon. 

Dalmont, Paris, 1856. [40 

F. Dullien. Porous Media: Fluid Transport and Pore 
Structure. Academic Press, San Diego, 2nd edition, 1992. 

C. Mei and J.-L. Aurialt. The effect of weak inertia on 

flow through a porous medium. J. Fluid Mech., 222:647- [41 
663, 1991. 

G. Massarani. Problemas em sistemas particulados. Ed. 
Edgard Blucher Ltda, Ri'o de Janeiro, 1984. 

P. Forchheimer. Wasserbewegung durch Boden. Zeit. [42 
Ver. Deut. Ing., 45:1781-1788, 1901. 

J. Andrade, U. Costa, M. Almeida, H. Makse, and 

H. Stanley. Inertial effects on fluid flow through disor- 
dered porous media. Phys. Rev. Lett., 82(26) :5249-5252, [43 
1999. 

E. Skjetne and J. Auriault. High-velocity laminar and 
turbulent flow in porous media. Transp. Porous Media, [44 
36:131-147, 1999. 

J.-L. Auriault, C. Geindreau, and L. Orgeas. Upscal- [45 
ing Forchheimer law. Transp. Porous Media, 70:213-229, 
2007. 

J.-C. Wodie and T. Levy. Corection non lineaire de la loi 
de Darcy. C.R. Acad. Set. Pans, 312:1573-161, 1991. [46 
M. Firdauss, J.-L. Guermond, and P. Le Quere. Nonlin- 
ear corrections to Darcy's law at low Reynolds numbers. 
J. Fluid Mech., 343:331-350, 1997. [47 
S. Chen, H. Chen, D. Martinez, and W. H. Matthaeus. 
Lattice Boltzmann model for simulation of magnetohy- 
drodynamics. Phys. Rev. Lett., 67(27) :3776-3779, 1991. [ 
P. L. Bhatnagar, E. P. Gross, and M. Krook. Model for 
collision processes in gases. I. small amplitude processes 
in charged and neutral one-component systems. Phys. [49 
Rev., 94(3):511-525, 1954. 

J. Sterling and S. Chen. Stability analysis of lattice Boltz- 
mann methods. J. Comp. Phys., 123(1): 196-206, 1996. [50 
Y. H. Qian. Fractional propagation and the elimination 
of staggered invariants in lattice-BGK models. Int. J. 
Mod. Phys. C, 8(4):753-761, 1997. [51 
M. C. Sukop and D. T. Thorne (Jr.). Lattice Boltzmann 
Modeling, An Introduction for Geoscientists and Engi- 
neers. Springer, second edition, 2007. 

D. d'Humieres, I. Ginzburg, M. Krafczyk, P. Lallemand, 
and L.-S. Luo. Multiple-relaxation-time lattice Boltz- 



mann models in three dimensions. Phil. Trans. R. Soc. 
Lond. A, 360(1792):437-451, 2002. 

P. Lallemand and L.-S. Luo. Theory of the lattice Boltz- 
mann method: dispersion, dissipation, isotropy, Galilean 
invariance, and stability. Phys. Rev. E, 61(6):6546-6562, 
2000. 

I. Ginzburg, F. Verhaeghe, and D. d'Humieres. 
Two-relaxation-time lattice Boltzmann scheme: about 
parametrization, velocity, pressure and mixed boundary 
conditions. Comm. Comp. Phys., 3(2):427-478, 2008. 
I. Ginzburg, F. Verhaeghe, and D. d'Humieres. Study of 
simple hydrodynamic solutions with the two-relaxation- 
times lattice Boltzmann scheme. Comm. Comp. Phys., 
3(3):519-581, 2008. 

S. Chapman and T. G. Cowling. The mathematical the- 
ory of non-uniform gases. Cambridge University Press, 
second edition, 1952. 

D. A. Wolf-Gladrow. Lattice-Gas Cellular Automata and 
lattice Boltzmann models. Springer, 2005. 
D. Kandhai, D.-E. Vidal, A. Hoekstra, H. Hoefsloot, 
P. Iedema, and P. Sloot. Lattice-Boltzmann and finite 
element simulations of fluid flow in a smrx static mixer 
reactor. Int. J. Numer. Meth. Fluids, 31:1019-1033, 1999. 
S. Srivastava, K. Yazdchi, and S. Luding. Meso-scale 
coupling of FEM/DEM for fluid-particle interactions. In 
preparation, 2012. 

H. Langtangen, K.-A. Mardal, and R. Winther. Numer- 
ical methods for incompressible viscous flow. Adv. in 
Water Res., 25(8-12):1125-1146, 2002. 
K. Yazdchi, S. Srivastava, and S. Luding. Microstructural 
effects on the permeability of periodic fibrous porous me- 
dia. Int. J. of Multiphase Flow, 37(8):956-966, 2011. 
Z. Guo and T. Zhao. Lattice Boltzmann model for in- 
compressible flows through porous media. Phys. Rev. E, 
66:036304, 2002. 

Q. Zou and X. He. On pressure and velocity boundary 
conditions for the lattice Boltzmann BGK model. Phys. 
Fluids., 9(6): 1591-1598, 1997. 

M. Hecht and J. Harting. Implementation of on- 
site velocity boundary conditions for D3Q19 lattice 
Boltzmann simulations. J. Stat. Mech.: Theor. Exp., 
2010(01):P01018, 2010. 



